%Matlab code to generate figures 1-3 in the manuscript "Putting the
%significance of spectral peaks on the level: implications for the
%1470-yr peak in Greenland δ18O" by P. Huybers, Feb 2022.  This
%file was run on Matlab release R2021b and requires the signal
%processing toolbox. 

clear;
addpath Toolbox;

%%load data
GISP2;
timeo=G(1:end,3)/1000;
Oo=G(:,2);
pl=find(timeo>=20 & timeo<=50);
O=Oo(pl);
time=timeo(pl);
mdt=mean(diff(time));
pl=find(timeo>=0 & timeo<=70);
Oo=Oo(pl);
timeo=timeo(pl);

%interpolate to uniform intervals
ti=time(1):mdt/10:time(end);
Oi=interp1(time,O,ti); 
Oi=smoothPH(Oi,11); 
time=ti(1:10:end); 
O=Oi(1:10:end);
dt=diff(time(1:2));

%%spectral analysis 
opts.p=1;
opts.q=0;
opts.nw=2;
[P1 s coeffs]=specw(detrend(O),dt,opts.nw,opts);

opts.p=1;
opts.q=1;
[P2 s coeffs]=specw(detrend(O),dt,opts.nw,opts);

opts.p=3;
opts.q=2;
[P3 s coeffs]=specw(detrend(O),dt,opts.nw,opts);

vP(1)=var(P1.prewhiten);
vP(2)=var(P2.prewhiten);
vP(3)=var(P3.prewhiten);

%F after removing the 1/(1470 yr) peak
pl=find(s<0.8/1.470 | s>1.2/1.470);
vP2(1)=var(P1.prewhiten(pl));
vP2(2)=var(P2.prewhiten(pl));
vP2(3)=var(P3.prewhiten(pl));

dof=2*(2*opts.nw-1);
k=dof/2;
Vexp=psi(1,k);  %expected variance 
F=vP./Vexp;

so=1/1.470;
ds=0.1/1.470;
[dummy plso]=min(abs(s-so));
mag1470(1)=P1.prewhiten(plso); 
mag1470(2)=P2.prewhiten(plso); 
mag1470(3)=P3.prewhiten(plso); 
save EX_GISP.mat mag1470 s plso vP;

%%plot time series
figure(1); clf;
subplot(5,1,[1 2]); cla; hold on;
plot(timeo,Oo,'color',[0.6 0.6 0.6],'linewidth',1); axis tight;
plot(time,O,'k','linewidth',1); axis tight;
set(gca,'xdir','reverse');
font(gca,14); 
h=xlabel('time (ky BP)'); font(h,14);
h=ylabel('\delta^{18}O'); font(h,14);
h4=axis; dhx=diff(h4(1:2)); dhy=diff(h4(3:4));

%%un-prewhitened spectral estimate
figure(2); clf; 
subplot(3,3,[1 2 4 5]); cla; hold on;
xline(so,'k--','linewidth',1);
pls=1:length(s);
plot(s(pls),P1.initial(pls),'k','linewidth',1); 
plot(s(pls),P1.pre_fit(pls),'r','linewi',1); 
plot(s(pls),P2.pre_fit(pls),'b','linewi',1); 
plot(s(pls),P3.pre_fit(pls),'color',[0.3 0.7 0.3],'linewi',1); 
set(gca,'xscale','log'); 
axis tight;
h=xlabel('frequency (1/ky)'); font(h,14); 
h=ylabel('power spectral density (log units^2/\Delta s)'); font(h,14); 
font(gca,14);
plot(3*[1 1],[1 1+coeffs.conf],'k','linewidth',1.5); 
plot([2.85 3.15],[1 1],'k','linewidth',1.5); 
plot(4*[1 1],[1 1+coeffs.confb],'k','linewidth',1.5); 
plot([3.8 4.2],[1 1],'k','linewidth',1.5); 

mx=P1.initial;
pls=find(s>=so-ds & s<=so+ds);
fp(1)=sum(exp(mx(pls)))/sum(exp(mx));
mx=P1.prewhiten;
fp(2)=sum(exp(mx(pls)))/sum(exp(mx));

figure(3); clf;
for rep=1:3,
    if rep==1, P=P1; col='r'; lab='ab'; end;
    if rep==2, P=P2; col='b'; lab='cd'; end;
    if rep==3, P=P3; col=[0.2 0.7 0.2]; lab='ef'; end;
    
    subplot(3,3,[1 2]+(rep-1)*3); cla; hold on;
    xline(so-ds,'k--','linewidth',1);
    xline(so+ds,'k--','linewidth',1);
    plot(s,P.prewhiten,'k','linewidth',1);
    plot(s,zeros(size(s)),'color',col,'linewi',1); 
    plot(s,coeffs.conf+zeros(size(s)),'--','color',col,'linewi',1); 
    plot(s,coeffs.confb+zeros(size(s)),'-.','color',col,'linewi',1); 
    axis tight; h4=axis; axis([h4(1:2) -3 2.5]); 
    h=xlabel('frequency (1/ky)'); font(h,13); 
    h=ylabel('leveled density (\psi)'); font(h,13); 
    font(gca,13);
    h4=axis; dhx=diff(h4(1:2)); dhy=diff(h4(3:4));
    h=text(h4(1)+dhx*0.01,h4(4)-dhy*0.03,lab(1)); font(h,14);

    subplot(3,3,rep*3); cla; hold on;
    xi=-4:0.01:3;
    Fk=ksdensity(P.prewhiten,xi);
    plot(xi,Fk,'k','linewidth',1);
    k=coeffs.dof/2;
    plot(xi,loggamma(xi,k,1/k),'color',col,'linewidth',1)
    h=ylabel('distribution'); font(h,13); 
    h=xlabel('leveled density (\psi)'); font(h,13); 
    font(gca,13);
    axis tight;
    h4=axis; dhx=diff(h4(1:2)); dhy=diff(h4(3:4));
    h=text(h4(1)+dhx*0.01,h4(4)-dhy*0.03,lab(2)); font(h,14);
    
    [rep mean(P.prewhiten>coeffs.conf), mean(P.prewhiten>coeffs.confb)];

    
end;

return;
for nex=1:3,
    figure(nex);
    pause(0.5);
    eval(['print -depsc fig',num2str(nex),'.eps']);
    eval(['!pstopdf fig',num2str(nex),'.eps']);
    eval(['!open fig',num2str(nex),'.pdf']);
end;


